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Abstract 

We consider the relativistic quantum mechanics of a two interacting fermions system. We 
first present a covariant formulation of the kinematics of the problem and give a short outline 
of the classical results. We then quantize the system with a general interaction potential and 
deduce the explicit equations in a spherical basis. The case of the Coulomb interaction is 
studied in detail by numerical methods, solving the eigenvalue problem for j = 0, j = 1, j = 2 
and determining the spectral curves for a varying ratio of the mass of the interacting particles. 
Details of the computations, together with a perturbative approach in the mass ratio and an 
extended description of the ground states of the Para- and Orthopositronium are given in 
Appendix. 
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"nI" . 1 Introduction 

p 

The relativistic two body problem and the properties of bound states in relativistic quantum 
mechanics have a very long history and a really large number of solutions have been devised. 
Some of these proposals stress the properties of covariance of the model, others assume as a 
starting point the constraint dynamics or arise directly from a classical symplectic context, others 
determine effective interactions from field theory in order to get results to be compared with 
phenomcnological data. The bound state spectra are considered with particular interest, being 
rather difficult to treat in a pure quantum field theory framework. A historical overview of the 
subject can be found in pQ. The research, however, has actively continued up to present days and 
several papers keep being published both on the mathematical and the phenomenological aspects 
of the problem (see e.g. and references therein). 

Some years ago the present authors too gave a contribution to this discussion, proposing a 
Lorentz invariant classical model for the two body Coulomb and gravitational interactions where 
the relativistic covariance is nonlinearly realized |14j . The paper was written in a plain coordinate 
language but it was rather strongly influenced by geometrical and symplectic ideas. These were 
developed elsewhere from a mathematical viewpoint and in a more general context, using group 
theoretical ideas for the treatment of classical spinning particles, [TJ], or studying rigorously the 
mathematics of skew-commutative phase spaces, |16| . We could however observe that a peculiar 
property of the coordinate system determined in Q3i ~ and briefly recalled in the following Section 
2 for the sake of completeness - is that the relative time coordinate disappears from the dynamics, 
and a symplectic reduction of the phase space is made possible, |17j : the relative time plays the 
role of a gauge and therefore it can be given an arbitrary form. The procedure permits then, if 
needed, to reconstruct the two separate world-lines in the Minkowski space, whose structure and 
correlation are obviously dependent upon the chosen form of the relative time gauge function. 

The present paper makes a revisitation of our previous model and deals with its canonical 
quantization, assuming that the interacting particles are fermions with arbitrary masses. We 
therefore develop our arguments and we find the wave equation for the two fermion system in a 
framework of pure relativistic quantum mechanics, without any reference to second quantization 
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techniques: this is, probably, the simplest way to transfer the classical results to the quantum 
context and to obtain a complete control of the symmetry and covariance properties of the theory. 
The wave equation is determined in the presence of a general scalar potential, but the explicit 
results are focused on the Coulomb interaction, with the purpose of calculating the spectrum of 
the lowest bound states and in particular of determining the dependence of the levels upon the 
mass ratio of the interacting particles for any value of the total angular momentum. It proves 
that the equations and the results are continuous in the mass ratio, so that they provide a unified 
scheme ranging from the Positronium (two equal masses) , to the Hydrogen atom and finally to the 
Dirac limit (one of the masses becomes infinite). Although written in a rather compact form, the 
equations we find are not easy to study analytically. This means that no simple formula for the 
energy levels is available as in the case of the Dirac equation and numerical procedures are therefore 
necessary if we want to determine the spectroscopy of the system. Appropriate plots will illustrate 
the results obtained: from these it will be possible to establish a regular and simple behavior of 
the energy levels vs. an appropriate mass ratio parameter. We find, for instance, that for two 
finite masses the degeneracy of the Dirac levels 2s 1 / 2 -2p 1 / 2 is removed and for the Hydrogen 
atom we reproduce a small fraction of the entire Lamb shift: it is interesting to observe that the 
correct kinematics gives a not negligible fraction of the splitting. The continuous dependence of 
the spectrum upon the ratio of the component masses suggests also in a natural way a perturbative 
treatment in the neighborhood of the Dirac limit. Except for the Positronium, the main systems of 
physical interest - as the Hydrogen atom and, at a lesser degree, the Muonium and the mesic atom 
- have a very small mass ratio and can be treated accordingly. We have shown how to develop 
this framework and we have applied it to the ground state of the Hydrogen atom. 

A comparison with the results of ref. is in order. Indeed, the equations found in our 

quantum mechanical framework agree with those deduced by different methods in the quoted ref- 
erences, although a bit of work is usually required to reduce the latter. Concerning the results, we 
may observe that all those papers develop a perturbative approach in the fine structure constant 
up to or; complete numerical results are given for equal masses, jH], and for different masses with 
a vanishing total angular momentum, |l()j . In |Hj we find a perturbative formula, corrected up to 
a 4 , expressing the Lamb shift due to the Coulomb interaction. The results of that formula are 
practically coincident with our numerical results for the Hydrogen atom, but appear to be less 
accurate when the masses of the two fermions become closer and closer: the comparison suggests 
that, at a high degree of accuracy, the perturbation expansions must be treated with some caution 
and numerical computations appear to be necessary and more reliable. Complete numerical results 
for equal masses and values zero and one of the angular momentum, are found in [jj]. The computa- 
tions are carried out by finite element methods and they are in complete agreement with our results 
obtained by the numerical solution of a singular boundary value problem. The authors observe 
that the Positronium presents an inversion in the energy of the two levels indicated as 2sxy 2 and 
2p 3 / 2 in the Dirac spectroscopy, the first being higher than the second. As a matter of fact we have 
drawn the complete plot of these terms for any mass ratio, showing their actual crossing. Finally, 
in the two papers of reference |10|. quantum field theoretical methods are used to determine the 
equations for fermion systems under the Coulomb and the complete electromagnetic interactions 
respectively. The equations are derived by defining an appropriate non conventional vacuum and 
by studying the conditions to be imposed on the coefficient functions of the states with a fixed 
number of particles in order to be eigenstates of the particle - electromagnetic field Hamiltonian. 
The interest is especially addressed to the properties of the system related to the variation of the 
coupling constant a and the perturbative aspects in terms of a are widely developed. Complete 
numerical results are given for a vanishing value of the angular momentum and perturbative results 
up to a 4 for non vanishing angular momentum. The wave equation we derive for the two fermions 
agrees with the results of upon restriction to the static Coulomb interaction. 

The two fermion relativistic bound states are relevant not only for electrodynamical purposes 
and atomic physics, but also in the framework of phcnomcnological meson spectroscopy. The 
limited progress of the QCD lattice program for the mesons has stimulated the interest in new 
treatments of the two body quantum systems that could replace the non relativistic Quark Model 
1 1 1 j and predict the masses of the mesons, from lighter to the heavier ones, assuming from the 
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experimental data the least number of values for physical quark parameters. Numerical approaches 
have devised new methods of solving the Bethe-Salpeter equation with interactions derived from 
QCD effective potentials, see e.g. |12|. New wave equations from the quantization of constraint 
theories have been proposed in a long series of papers, see [5], and references therein. The 
authors start from two coupled Dirac equations; using the Poincare symmetry and some super- 
symmetric ideas, they are able to determine the general structure of the interactions compatible 
with the constraint scheme. The application to the phenomenology of mesons gives a very good 
agreement with the experimental data. 

We have written this paper trying to make it self-consistent. Some details of the computations 
have been explained in the Appendixes. A brief summary is as follows. In Scction|2]we recall the 
basic definitions and we establish the kinematics and the dynamics of the system together with 
their invariance properties and their validity in any reference frame. In Section|2|the quantization 
procedure is explained and applied to the relativistic two fermion problem. In Section 0] we deduce 
the equations for an arbitrary scalar potential. Given the general form of the state vector with 
definite energy, angular momentum and parity - whose explicit expression involves eight unknown 
functions and is reported in Appendix A - we give a short, transparent and comprehensive de- 
duction of the equations for the unknown functions and we specify them to the particular case of 
the Coulomb interaction. The computations leading to the final result are sketched in Appendix 
B. The Dirac limit is treated in Section [SJ here we present a procedure that permits to recover 
continuously the Dirac equation for the finite mass constituent when the mass of the other goes to 
infinity. Finally, in Section [H] we discuss the numerical results with j = 0, j = 1 and j = 2 for the 
whole range of masses: the results for j = should be compared with a first order perturbative 
approach in the ratio of the constituent masses, described in Appendix C. Some few details of the 
numerical procedure are given in Appendix D. In the final Appendix E we report on some details 
of the ground state of the Positronium. In particular we calculate the wave function of the Para- 
positronium and Orthopositronium, we show that the non relativistic limit correctly reproduces 
the Schrodinger equation for the Coulomb problem and we make some comments on the hyperfine 
splitting. 

In the paper we shall use a unit system with h = c = 1, and for the presentation of the results 
we also assume the electron mass equal to unity. 



2 The classical problem and the dynamical variables 

We first outline the procedure for two classical scalar particles with given masses mi, i = 1 , 2, whose 
individual Minkowski coordinates are denoted by and whose canonical conjugate momenta are 
. Starting from the collective coordinates 
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it has been shown in [T2|, that a good set of canonical conjugate variables is provided by the 
following definitions: 
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where rf iV is the Lorentz metric tensor, e a b c is the three dimensional skew-symmetric tensor and 
where the sum over repeated indices must be assumed. In equation l|2.2|) the tensor (^t, a — 0, 3) 
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given by 



aV ' la VP 2 [P a + VP 2 } ' 

eft(P) = P^/Vp 2 (2.3) 

realizes the Lorentz transformation to the P a — reference frame and therefore it satisfies the 
identities 

V^^(P)ef 3 (P) = T] a p, 

r, a psZ(P)e» (P) =rT. (2-4) 

Indeed it can be observed that both r a and q a are Wigner vectors of spin one, as well as Z a has the 
structure of a position vector of the Newton- Wigner type for a particle with angular momentum 
L a , where 

L a = £ab c r b q c . (2.5) 

We can remark that the variables l|2.2J) arise quite naturally as global and relative coordinates 
in a two body Poincare invariant dynamics constructed by using the algebraic and coalgebraic 
properties of the Weyl homogeneous spaces. Indeed a four component position operator for each 
constituent is built in terms of the generators of the corresponding Weyl algebra, whose non trivial 
cohomology permits to deduce global and relative operators. The breaking of the scale invariance 
leaves the resulting dynamical system Poincare symmetric, |15j . This procedure, although different 
in the algebraic assumptions and in the results, presents clear analogies with the classic paper |18| . 
(see also subsequent literature, e.g. ^H]), where the two body theory was for the first time 
coherently based on the algebra of the Poincare generators corresponding to each particle. 

Using (|2.2I) , the mass shell conditions for each of the two particles read 

P(0 = i p2 + H i+1 + ~ = m i ■ ( 2 - 6 ) 

from which 

VP 1 q= \{m\-m\), 

\P 2 + 2q 2 - 2q a q a = m 2 + m 2 . (2.7) 



Therefore the total mass A = yP 2 results in 

y q a q a +m 2 \ +(q a q a + ml\ =A, (2.8) 

while the variable q can be fixed, generating a symplectic reduction of the phase space. Its conjugate 
coordinate r - the relative time coordinate - is cyclic and assumes the character of a gauge function 
that is chosen a posteriori in order to recover the complete Minkowski description for the two 
particles. In particular cases it could be useful to fix r = 0, but in principle there is no necessity 
of requiring such condition. 

It is finally straightforward to construct a Lorentz covariant dynamics by introducing nontrivial 
interaction terms depending only upon q a and r a . One of the simplest choices is to add scalar 
potentials that are functions of the Lorentz scalar r = (r a r a y^ 2 - Hence, a relativistic two-body 
system interacting by means of a potential V(r) is described by 

yqaqa + mlj +{q a q a +ml\ =h(r). (2.9) 

with 

h(r)=X-V(r) (2.10) 
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The equation of the orbits for the Coulomb interaction, V(r) — —a/r, is obtained from the 
canonical equations deduced from l|2.9|) . Introducing the variables u and defined by 



u = A . 

r 



cos = — 



(2.11) 



and denoting by L the value of the conserved angular momentum, the equation reads 



cW 
du 



(AL 2 /a 2 )u 2 (X-u) 2 -2(m 2 + m 2 )u 2 + (m 2 -m 2 ) 2 {Atf/a 2 ) 1 ' 2 (2.12) 



-1/2 



and it can be integrated in terms of elliptic functions (see ^2] for details). In the case of equal 
masses, mi = m 2 — to, the solution for the orbits can be given in terms of elementary functions. 
For instance, for a 2 < AL 2 - the most relevant case from our point of view as it gives bounded 
orbits - the solution is 



4L 2 - a 2 



4L 2 A 2 + 4 (a 2 - 4L 2 ) to 2 ) cosh [ 1 - a 2 /AL 2 } 



211/2 



A a 



(rara) 1 / 2 

and corresponds to the motion of a particle in an external Coulomb held (see e.g. |19| 



(2.13) 



3 The quantization of the free two fermion system 

The next step is the quantization of equation (|2.9|l with V(r) = 0. The scalar quantization of 
the double square root Hamiltonian has received many attentions and has been discussed also 
recently in terms of functional inequalities, |21|. However, since most of the physically interesting 
situations involve fermionic components, we assume particles of that nature and we quantize the 
system accordingly. 

The Dirac operators corresponding to each single particle are 

Di = (^-P/i + cfc) 7(1) - Til , D 2 = (^P„ - 9 M ) 7(2) - m 2 (3.1) 
where we have introduced the following tensor products of gamma matrices 

Tfo = 7 M ® 14 , 7f 2) = 14 ® r , (3.2) 

where I4 is the unity matrix in four dimensions. 

The operators Di and D2 can be rewritten in terms of the canonical set l|2.2|l . At the same time 
the constant 7^-matrices will be recast in terms of the new set of matrices {*y(P), ("f a (P))a=i,3}, 
where 

7 (F) - £ £(P) 7 „ , la{P) = sZiPhu (3.3) 
The two Dirac equations (|3.1|l become then 

~ A 7(i) + 97(1) - q a l(i) a = m 1 , ^A7 ( 2) -97(2) +9«7(2) a =m 2 . (3.4) 

We observe that the square of 7(j) gives the unity matrix. We next multiply the equations l|3.4|l 
on the left by 7m and 7(2) respectively and we obtain 

2 A + g-tfo 7(i)7(i) „ = 7(1) "M > ^ A - q + q a \ 2 ) 7( 2 ) Q = 7(2)TO 2 . (3.5) 
By summing and subtracting we finally get 

A = 3o(7(i) 7(i) - 7(2) 7(2) ) + 7(1)^1 + 7(2)^2 (3.6a) 
9= 5 80(7(1) 7(i) a + 7(2) 7(2) ) +7(1) mi -7(2) "12 (3.6b) 
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The right hand sides of the two equations (|3.6[) are commuting. Multiplying those equations term 
by term, we get 



\ " 2 2 
\q = m 1 — m 2 



(3.7) 



and we see that the variable q remains fixed, in complete agreement with the classical symplectic 
reduction. In fact the canonical conjugate relative time coordinate r is cyclic again: this is signified, 
in particular, by the Lorentz scalar identity for the phase of plane waves 



q a r a 



(3.8) 



i ft \~i \~/fi 

We finally observe that the definition (|3.3|) is actually a Lorentz transformation on the 7^ four 
vector determining a unitary transformation. Therefore, as long as P is conserved, the matrices 
(|3.3() will be represented by the usual 7 matrices. The different behavior under Lorentz transfor- 
mations is however signified by the use of the notation Ij3.3|l . In conclusion equation Ij3.6a|l is the 
Lorentz-invariant equation for the two-fermion system. It is also straightforward to calculate the 
explicit expressions for the sixteen eigenvalues A and the result is rather obvious. In fact we have 
four singular values 



A = ± I q a q a 



± 



1/2 



± (q a q a 



(q a q a 



1/2 



(<7a q a 



1/2 



1/2 



(3.9) 



each one with multiplicity four. Defining 

M = mi + m,2 , [i 



mi — mi 



P = P/M . 



(3.10) 



we make a linear transformation on the tensor product of the two spinor spaces that diagonalizes 
the system at rest, in such a way that the four singular values (|3.9|) appear in the order M, —A/, 
— ji and fi. We find it also convenient to perform a further linear transformation that, in each four 
dimensional eigenspace, diagonalizes the square and the third component of the total spin 



S = I 4 



(3.11) 



<r being the Dirac spin. For each singular value of the mass, the diagonalizcd spin will then be 
ordered with the triplet always following the singlet. 

The wave equation will now be obtained by applying the canonical quantization rules to the 
equation (|3.6a|l . In view of the bases we have chosen, the natural coordinates to be used are the 
spherical ones, in terms of which the free Hamiltonian operator reads 



#0 = 

where each matrix element is actually an ! 

I 4 \ 



Jm = M 



( Jm T~(-o \ 
\ Ho Jft J 
x 8 block. Explicitly 
X I 4 
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J/t = M 
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Xo 


X- 
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x+ 


Xo 


-x+ 
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-Xo ) 



(3.14) 



G 



The matrix elements in <|3.14[1 are the spherical operators 

X± = -2-^ 2 {±q x +tq y ) 1 X = q z , (3.15) 

where q a — > —id/dx a ■ 

We finally remind that the global parity transformation is given by the product of orbital and 
internal parity transformations. In our picture the internal parity is given by 




(3.16) 



It is straightforward to verify that the global angular momentum J = L + S and the parity 
are conserved. Together with A they provide a classification of the states of the global symmetry. 
Concerning the problem of canonical realization, Poincare representations and invariant scalar 
products we refer to the paper (201 ■ ^ is also straightforward to observe that the conservations 
properties keep holding in the presence of interactions depending upon the Lorentz scalars r a r a , 
q a q a and q a r a - this fact will be used in the next section to determine a reduced set of differential 
equations and to state the boundary value problem in the space of the square-integrable functions 
on the relative R 3 . 



4 The spectral problem for an interacting system 

The construction of the states with assigned angular momentum (j,m) and given parity (— ) J or 
(— y +1 is done as usual, by multiplying the contributions coming from the composition of orbital 
and intrinsic angular momenta by functions of r. The expressions of the states and \f r _ are 
reported in equations (|A.1|I and (|A.3(1 given in Appendix A. 

The eigenvalue problem comes out when we try to determine the eight unknown function 

Oi(r), bi(r), Cj(r), d;(r), (* = 0,1) (4.1) 

in the state \& by solving the homogeneous equation 

(fl o -/i(r))* = J (4.2) 

with h(r) as in (|2.1()(l . Since the parity and the angular momentum are conserved, equation l|4.2J) 
produces two different spectral problems to be discussed separately. 



4.1 The spectral problem for the state \P + 

When substituting ^ = ^ + in (|4.2|l . we obtain a system of equations given by the vanishing of the 
coefficients of the different spherical harmonics in each component of the resulting vector. This 
system is composed of thirty-four first order differential equations, but, as one should expect, only 
eight of them are independent. The detailed expressions for the eight equations are written in 
Appendix B. 

It turns out that appropriate changes of the initial variables can help in giving the system a 
much simpler and readable form. To this purpose, we first define the sums and differences 

s±(r) = s (r) ± si(r) (s = a,b,c 1 d) (4.3) 

and we consider the following linear combinations of the c±(r) and d±(r): 

t \ Vjc+(r) - yg + ld+(r) 
u +\ r ) = . 

v 1 \/TT27 
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u_(r) = — 
v+{r) = - 
V-(r) = — 





VH 


-2j 


Vj- 


hlc+(r) 


+ Vj^+W 




Vn 




Vj- 


h lc_(r) 


+ VJ^-W 




VI 4 


-2j 



(4.4) 



A rather simple formulation for the system can then be given in terms of a± (r) , b± (r) , u± (r) and 
v± (r) . Indeed, from a straightforward computation, we find that the eight independent equations 
split into four algebraic relations and four first order differential equations. The former read 



Vj(j + l)a+(r) - - (h(r) v+(r) +[*V-(r)) = 



Vj'C?' + 1) b-(r) + - 0*«+(r) + h(r) u_(r)) = 



Ma + (r) -h(r) a_(r) =0 
/i(r) 6+(r) - ML(r) = 
while the differential equations are 



(4.5) 



<Mr) , V^M Mr)+ (zM!M Mr) = 



c/r 



r /i(r) 



2ft(r) 



7 .7 



d6_(r) jh_(r)\ vJTT 



dr r 
du+(r) M-h(r) 
~~dr 2 

dv-(r) M — h(r) 



a+(r) + a-(r) ) + - u + (r) — - — - + v + (r) = 
r r 



b+(r)+b(r) ) - vSEES tt-(r) + i«_(r) = 
r r 



(4.6) 



dr 2 

By means of (|4.5f) we eliminate the four variables that are not differentiated, obtaining a 
system of four first order differential equations in the four unknown functions a+(r), b—(r), u+(r) 
and V-(r). If we arrange the latter respectively as the four components of a vector Y(r) = 
t (yi{r),y2(r),y3{r),y4(r)), the system is given the compact form 

dY(r) 



dr 



+ MY(r) =0, 



(4.7) 



where M is a matrix with general structure 



M = 






E(r) 


F(r) 





E(r) 


1 

r 





-F(r 


G(r) 





2 

r 


E(r 





-G(r) 


£(r) 


1 



(4.8) 



The explicit expressions of E(r), F(r) and G(r) for the even case are the following ones 



E c (r) = 



Vi(j + I)M 
r h(r) 



G c (r) = k M , , 



/i 2 (r 

r 2 Af 2 + 4j(j + I) 
r 2 /i 2 (r) 



(4.9a) 
(4.9b) 
(4.9c) 
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Specifying the matrix elements (|4.9|) to the case of the Coulomb interaction, where 



h(r) = X + a/r, (4.10) 
and where a is the fine structure constant, we have 



E c ,cou\(r) = ■ (4.11a) 

Ar + a 

W) = - (Ar 2 + rf:f 2 ( 4 - Ub ) 
2r (Ar + a) 

(Ar + a) 2 -Af 2 r 2 -4j(j + l) 

G e ,Coui(r- = 7T - s 4.11c 

2r (Ar + a) 

Some remarks are in order. First of all, from (|4.8|) and 1)4. 9|) we see that the transformation 

A* ~ *• — Mi 2/iM -» -2/iM i 2/3 W ->■ -2/3(r) (4.12) 

is a symmetry of the system, corresponding to the fact that the ordering of the two particle is 
irrelevant. 

In the second place, from l)4.9a)) . we see that E e (r) — either for j = or /i = 0. In these 
cases the fourth order system splits into separate second order subsystems, one for the functions 
j/i (r) and 2/3 (r), the other for 2/2 (r) and 2/4 (r). The subsystems, in turn, can be presented as the 
following two independent second order differential equations 



' yiW + (J ~ ^FW~) ^ 2/1 W " F W GW yi(r) = ° (4 - 13a) 



d 2 | /2 ^ F ( r )\ d 

<ir 2 \r -FX*") 

d 2 /2 T^K d \ 

^ + (- r - %r) * »w ^ ( w + G(r) ) m{r) = (413b) 

As discussed in Sectional the Dirac limit shows that for j = only (|4.13a|l makes sense. 

Although the linear system 1)4.7)1 and, in particular, the equation (|4. 13a|) do not look extremely 
complicated, they still are difficult enough to prevent the possibility of an easy analytical solution. 
For instance, according to a somewhat old fashioned - but actually still the most up to date - 
classification of differential equations according to the number of their elementary singularities (for 
definitions and detail we refer to the classical text [22]), equations 1)4.13)1 with the Coulomb coeffi- 
cients (|4.11|l present eight elementary singularities, while analytical properties of the solutions have 
been studied and classified for differential equations having at most six elementary singularities. 
Numerical techniques are therefore necessary if we want to determine the spectral properties, both 
for equal and different masses. 

4.2 The spectral problem for the state 

The differential equations generated by the action of the Hamiltonian on the odd state are easily 
derived by observing that the transformation (|A.3[) amounts to changing mi into —mi, which 
means /i — * — M and M — > — jj,. The general structure for the matrix of the odd linear system, 
therefore, remains the one given in 1)4.8(1 . The matrix elements, on the other hand, are as follows: 
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and also the Coulomb explicit expressions are obtained from (|4.11|) by the exchanges (M, p) — > 
(— fi, — A/). The system decouples for j = 0, but, in contrast to what occurs in the even case, this 
circumstance is not realized when the particles have equal masses. 



5 The Dirac limit 

In this section we prove that (|4.7I4.8|) has the correct Dirac limit. This means that when the mass 
of one of the two interacting particles, say mi, tends to infinity, the system reduces to the Dirac 
equation for a charged fermion in an external potential. 

The limit is defined by assuming that the mass 777,2 = m e remains finite and gives the mass 
scale. The rescaling of the mass parameters of the system is then 

2pm e 2m e 

p = , M = , (5.1) 

1-p 1 - p 

We also take a dimensionless independent variable 

x = m e r (5-2) 
and define a coefficient 77 by means of the relation 

h(x) — (ri - v(x) + j + P ^ m e (5.3) 

where v(x) is the potential expressed in terms of x and divided by m e . It proves it convenient to 
define the shorthand 

rj(x) = r) — v(x) . (5-4) 

We now take the limit p — > 1 on the matrix elements. For the even system in the dimensionless 
variable x, we find 



For the odd system: 



E c {x) -» y/M + l)/x 

F c (x) -» -(l + v (x)) 

G e (x) -> n(x) - 1. (5.5) 

Eo(r) -» - VJ{j + l)/x 
F a (r) - (1 - r,(x)) 

G {r) -» (l + 7j(ar)). (5.6) 

Let us prove that the equations obtained by the use of l|5.5|l and (|5.(i|l are actually equivalent to 
pairs of Dirac equations |23| . that we write in our rescaled variables as 



(5.7) 



It is not difficult to realize that a mixing is necessary in order to decouple the fourth order system in 
two second order subsystems to be compared to the Dirac equation. This is physically clear, since 
from a description focusing on the angular momentum properties of the bosonic global system, we 
want to determine the equation for a fermionic object that constitutes only a part of the compound. 
The required linear transformation is generated by the orthogonal constant matrix 
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u 



VW+l 



VT+T 
Vj 






-Vj 



Vj VT+i \ 



Vj 



vm 



T 



-Vj 
o 
o 



(5.8) 



By defining 



Zi(x) 



4 

£ 



(U xyj(x) 



(5.9) 



it is straightforward to see that both in the even and the odd case the equations for the variables 
Z\(x) and 2:4(2) decouple from the equations for 2:2(2) and 2:3(2). The comparison with l|5.7[) 
shows that, the subsystems obtained in the even case are respectively the Dirac equations with 
k = — (J + 1) and k = j (here we must change 23(2) into —2:3(2)): these are the two allowed values 
of k for the states whose orbital momentum is j in the Dirac spectroscopy, 2.'! . In the odd case 
the decoupled systems agree respectively with the Dirac equations with k = j + 1 (here changing 
24(2) into —24(2)) and n = — j: these values correspond to orbital angular momenta equal to j + 1 
and j — 1. 

The previous discussion clearly shows that the complete treatment is continuous in the param- 
eter 1 — p : indeed, when this parameter is small, it can be effectively used computing perturbative 
corrections to the Dirac equation. Thus, although the general features of the spectrum - e.g. its 
boundedness from below - are clearly dependent upon the potential function ^(2) and can be stud- 
ied in each particular case by well established criteria 24 , the continuity in 1 — p of the differential 
operator implies the continuity in 1 — p of its spectrum that therefore shares the qualitative prop- 
erties of the limiting case. We report in Appendix C the calculation of the first order correction 
to the Hydrogen ground state and we find an excellent agreement with the numerical results. 



6 Discussion of the results 

In this section we present the results obtained. We display the tables containing the values of the 
variable p and the corresponding level with ten meaningful figures; the corresponding plots are 
then drawn. The energy parameter w that appears both in the tables and in the plots is related 
to the parameter r\ defined in (|5.4|l by the relation 

77= 1+-(1 + p) a 2 w (6.1) 

In terms of the eigenvalue A, this equation is equivalent to A = M + mna 2 w, where tur is the 
classical reduced mass rriR — mi raij(ra\ + 7712) . This means that the correction of energy levels 
due to the classical reduction to the center of mass frame is already accounted by the scale. Hence, 
the dependence upon p shown below is a genuine relativistic effect. 

The results will be given in atomic units, i.e. m e = 1 in addition to K = 1 and c = 1. The 
values of p with greater physical relevance are p = p e , corresponding to the hydrogen atom, p = p^, 
that refers to the mesic atom and p = 0, for the positronium. The value of w for p = 1 has been 
taken from the analytical formula for the Dirac spectrum. It is natural to begin from the ground 
state. We find that for any value of the two fermion masses the lowest energy level is degenerate 
with multiplicity two and the corresponding states are the first (i.e. with lowest eigenvalue) even 
state with j = and the first odd state with j = 1. We specify the results in the Table here below 
and we shall comment on the degeneracy later on. 
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p 


w + , j = 0, I 


P 


3=1, I 


u.u 


/looon^m no 


u.u 


zioooo^m no 
-.4yyyyouiuy 


0.2 


-.4999954766 


0.2 


-.4999954766 


0.4 


-.4999968737 


0.4 


-.4999968737 


0.6 


-.4999992025 


0.6 


-.4999992025 


Pn 


-.5000024182 


P» 


-.5000024182 


Pe 


-.5000066312 


Pe 


-.5000066312 


1.0 


-.5000066566 


1.0 


-.5000066566 



Here p e = 0.998911 corresponds to the hydrogen atom data and — 0.797576 refers to the mesic 
atom. The plot of the ground state vs. p is shown in Fig.l. The difference w + j = q, \ {Pe) ~~ 
u> +! j = o.i (1) = 0.254 • 10~ 7 should be compared with the value 0.253 • 10 -7 obtained in Appendix 
C by a perturbative calculation. 
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-0,500006 
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FIG. 1 The dependence of the degenerate ground state upon p : the levels w + j ■ — o, i and w — j = i i . 

The triplet given by the first even excited state with j = - corresponding to the state 2si/2 
in the Dirac limit - together with the even states that in the Dirac limit reduce to 2p 1 / 2 and 2p 3 / 2 
is presented in the first plot of Fig. 2. Their numerical values are respectively reported in the Table 
that follows: 



p 


3=0, II 


P 


W +, 3 = 1, I 


P 


W + i 3 = 1, II 


0.0 


-.1249996884 


0.0 


-.1250005200 


0.0 


-.1250002427 


0.2 


-.1249997840 


0.2 


-.1250006292 


0.2 


-.1250002030 


0.4 


-.1250000710 


0.4 


-.1250008727 


0.4 


-.1250001675 


0.6 


-.1250005493 


0.6 


-.1250012005 


0.6 


-.1250001864 


/V 


-.1250012098 




-.1250015984 


Pn 


-.1250002673 


Pe 


-.1250020750 


Pe 


-.1250020774 


Pe 


-.1250004151 


1.0 


-.1250020802 


1.0 


-.1250020802 


1.0 


-.1250004160 
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FIG. 2 The spectral curves for varying masses from the even and odd states. In the left plot (even states), circles 
represent the curve w+ ,j = o, n that reduces to the state 2s]y 2 in the Dirac limit, triangles the curve 
w+^j — i } i reducing to the state 2pj/2 and squares the curve Wj rt j = x, n reducing to 2p%/2- In the right 
plot (odd states), circles to refer to iu_ j = i i , triangles to ui _ ,j = o, I > squares to w — j=2 I ■ 

In the final Table we give the values for the three states that constitute the odd counterpart of 
the even triplet and that are shown in the second plot of Fig. 2. They arise from the three different 
values j — 0, 1, 2 of the angular momentum. 



p 


W ~, 3=0, I 


P 


= II 


P 


W ~, 3=2, I 


0.0 


-.1250007974 


0.0 


-.1249996884 


0.0 


-.1249999654 


0.2 


-.1250008487 


0.2 


-.1249997840 


0.2 


-.1249999834 


0.4 


-.1250010026 


0.4 


-.1250000710 


0.4 


-.1250000375 


0.6 


-.1250012592 


0.6 


-.1250005493 


0.6 


-.1250001276 


/V 


-.1250016134 


Pn 


-.1250012098 


Pn 


-.1250002520 


Pe 


-.1250020774 


Pe 


-.1250020750 


Pe 


-.1250004150 


1.0 


-.1250020802 


1.0 


-.1250020802 


1.0 


-.1250004160 



From the data we discover one more degeneracy: indeed the even state w + j = o. n an d the odd 
state w-j = i t ii, obviously equal in the Dirac limit, remain degenerate for any value of p. This 
occurrence, in the present case and for the fundamental state as well, is certainly not evident from 
a mathematical or computational point of view: indeed the equations that have been solved to 
find the even spectral curves and the systems that give the the odd ones are completely unrelated. 
We can however understand this apparently peculiar behavior by thinking of the non relativistic 
limit. Here, unless the spin-orbit coupling is explicitly included, the total spin and the orbital 
angular momentum are separately conserved. Moreover, if no spin-spin interaction is present, 
the relative orientations of the two spins between themselves and with respect to the orbital 
angular momentum are effective in determining the parity, but have no influence on energy of 
the state. The result is that the even and the odd states are degenerate with respect to parity. 
This degeneracy survives in the Dirac limit: although the spin-orbit interaction is now directly 
accounted for, one of the two constituent particles is infinitely heavy and thus at rest, so that only 
the spin of the light particle matters. The situation is different in our relativistic equation, when 
both constituents have finite mass and each of them gives its own contribution to the spin-orbit 
coupling: in fact we obtain different shifts for the even and odd terms, but for those states that in 
the non relativistic limit have a vanishing orbital angular momentum. In the particular case of the 
Positronium, we see that the two states iW-|-,j=o, i and W- t j — x t i correspond to Parapositronium 



13 



and Orthopositronium respectively. Their degeneracy could be removed by introducing the spin- 
spin interaction, responsible for the hyperfine splitting. We shall comment on this subject in 
Appendix E. 

Coming back to the figures, we observe that the solid curves are the parabolas obtained by the 
best fit of the data and they are in excellent agreement with the data themselves. The perturbative 
formula given in [8] for the Coulomb Lamb shift, 6w — (a 4 /12)(l — p)/(l + p) 2 , can reproduce 
a parabolic behavior only for p in the neighborhood of zero or unity. Only the second of these 
possibilities is actually realized: this shows the necessity of a more detailed numerical analysis 
when the masses become closer and closer. It is nevertheless interesting to observe the increasing 
impact of the relativistic effects with the decreasing of the mass difference. In particular it has to be 
noticed the crossing of the levels w j =o, n and W-\- t j = i, n as well as w _ j — 1. n and w - t j =2, i : 
this, for instance, reflects the necessity of the different method in use to classify the Positronium 
states with respect to the Dirac scheme. 

We conclude this discussion of the results with some observations on a possible formulation of a 
more realistic model involving a minimal coupling that accounts for the complete electrodynamical 
interaction. In fact this type of coupling can be arranged quite easily in our quantum mechanical 
framework and, if we study the case of a closed two-body system and we do want to keep a pure 
relativistic quantum mechanical description, we can try to express the vector potential itself in 
terms of the dynamical variables of the two particles. The Lienard-Wiechert potentials provide 
the most natural answer to the raised question. Unfortunately they introduce a delay in the 
propagation of the interaction that induces a considerable amount of difficulty in the solution of 
the corresponding wave equation and indeed the necessity of rendering its solution less tough is a 
major source of approximations, that, in most cases, produce not very satisfactory results (see [S] 
and references therein for a discussion on the subject). Indeed if one stops at the lowest order of 
approximation in the Lienard-Wiechert potentials, neglecting the the propagation delay, one finds 
an interaction energy 



that has to be substituted to V(r) in l)2.10[l . According to the Dirac prescription, when switching 
to quantum mechanics, the particle velocities Vu\ a will be substituted by the a^ a — — Ju)J(i) a 
matrices, so to obtain an invariant formulation of the Breit approximation for the electromagnetic 
interaction. The final wave equation reads 



The difficulties of a direct solution of equation l|6.3[l have been known for a long time (see |25 ] .|2S ] 1: 
a perturbative interpretation of the term involving the vector potential, averaged over the pure 
Coulomb wave function - very often at the non relativistic Schrodinger level of approximation -, is 
therefore what it is usually computed. The calculation can be found in |S] and produces a perfect 
cancellation of the a 4 -perturbative Coulomb splitting. From the previous discussion it is clear 
that the result should be tested numerically using the relativistic Coulomb wave functions for any 
value of the mass ratio. This numerical test is however beyond the purpose of the present paper. 
In any case - and as it was certainly to be expected - it appears that the acceleration and delay 
effects of the Lienard-Wiechert potentials - and therefore the radiative corrections - are essential 
for obtaining physical results. 



a a f 
V m (r) = - + - 
r 2r V 




(6.2) 





(6.3) 
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Appendixes 



A The state vectors 

We give here the explicit form of the 16 component state vectors of definite energy, angular mo- 
mentum (j,m) and even and odd parity with respect to the angular momentum, namely (— )■? and 
(— y +1 . We call them ^ + and ^_ respectively. The even state is then given by 

*+ = *(X M \ *(_- M \ *<"">, *M ) . (A.l) 

The four components is actually multiplets composed by a singlet and a triplet and will be indicated 



as: 



(A.2) 



where A = ±M , . According to a common practice and in order to deduce real differential 
equations in the following Appendix B, we shall introduce an imaginary unity in front of the 
coefficient functions Cj (r) and di (r) The explicit expressions of the state vector components are the 
following: 

^=Yl(0, 0) ao(r) 



(M) V] -m+ l^j + m j 

3^-1(0, 4>) M r ) 



/27VJTT 



AM) 



V27VJ+T 

Vto M) =^(0, 4>) 01 (r) 



0) 6o(r) 



v^v / T+7 



,y*(0, 0) &i(r) 



^(-M) = ^v^ ±T y , +i(M) bi(r) 
^=0 



/2j y2j^T 



YJn-iW, 0)*co(r) + -===—== y^.!^, ^)ido(r) 



V2j + 2V2j + 3 
V j - m + 1 yjj + m + 1 



Yl+\9, 4>)ido(r) 
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y m+1 (y,^ lCl (rj+ y m+lW 0j idl (rj 

where Y£(6, (j>) are the spherical harmonic functions. Concerning the state with opposite parity, 
it can be seen that the action of the parity transformation is simply the exchange of the first eight 
components with the second eight ones. Therefore, in order to determine the differential equations 
describing the odd states of the system, we shall use the state ^F- given by 




•-U °r + (A - 3) 

As observed in Section IV. B this amounts to changing the sign of the mass mi. 



B The original system 

As we said in Section (IV. A), in order to determine the system of differential equations describing 
the two body problem, we shall apply the Hamiltonian operator Hq, given in (|3.12|) . to the state 
vector IjA.ljl . In each component we require the vanishing of the coefficient of each different 
spherical harmonics. This leads to a system of thirty- four differential equations in the unknown 
functions cii(r) , bi{r) , Ci(r) and di(r) , (i = 0, 1) , that appear in the state vector. However, as 
expected, only eight of these differential equations are independent. They read 

+ J -^~) Mr) + 0l (r)) - V7 + T (± + £±i) (6 (r) - h(r)) 

+ ^2j + lc (r) (/x + /i(r)) = (B.la) 

(± + 1±1) ( ao(r ) + 0l ( P )) + vTTT (1 + (6„(r) - h(r)) 

\ar r / \ar r J 



v/2j + lci(r) (p-h(r)) = (B.lb) 



v/7+T - J) (a (r) + ai(r)) + Vi (- - (6 (r) - & x (r)) 



-V^7+ld (r)(/i + /i(r)) =0 (B.lc) 

V7 + T (J: - ~) (ao(r) + ar(r)) - V7 (J; - J) (&o(r) - &i(r)) 

+ V^7+ldi(r)(/i-/i(r)) = (B.ld) 

^ (i " S^) (co(r) + ci(r)) " v/7TT + ^) (do(r) + di(r)) 

+ y/2j + lao(r)(M-h(r))=0 (B.le) 
(J: - ^) (co(r) + Cl (r)) - V7TT (1 + i±^) (do(r) + di(r)) 
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-v/2j + lai(r) (M + h(r))=0 (B.lf) 

(c ° (r) - ci(r)) + ^ a + 3j ¥) (d ° (r) - di{r)) 



-y/2j + lbo{r) (M-h(r))=0 (B.lg) 

V7+1 (J: - ^) (co(r) - c x (r)) + VJ ( J: + ^) (do(r) - di(r)) 

- V57+16i(r) (M + ft(r)) =0 (B.lh) 

Once the definitions 1(4. 3() and 14.4fl are used, a straightforward, although somewhat lengthy 
computation, leads directly to 14.5(1 and 14.6(1 . 

C A perturbative calculation in the mass ratio 

We now present a perturbative approach to the ground level of the Hydrogen atom, taking as 
perturbation parameter the quantity 

e = 1 - p H = 1 - mp - me = 0.0010886411 (C.l) 

mp + m e 

As we have said in Section 0] for j = the system splits in two subsystems obtained by (|4.7I4.8[I . 
or, equivalently, by the second order equation ((4.13a(l . For later convenience in the development of 
the perturbative calculations we introduce the modified unknown functions <fi(r), 9(r) defined by 

yi(r) = ^5 0(r) , aa(r) - 6 {r) . (C.2) 

r r 

In particular we see that second order equation 14.13a(l written in terms of 4>(r) does not present 
the first derivative term: 



h^^l^W) ^)^)-|^-)^(r)=0 (C.3) 



For F(r) and G(r) we shall then assume their Coulomb explicit form 1(4. 9|) . Finally, the parameters 
are rescaled according to 15.115.3(1 and, correspondingly, we write the series expansions of the 
eigenvalue and the eigenfunction components in terms of e : 

A = ^Te"A n , y j {x)=Y J z n yf\x), (7 = 1,3). (C.4) 

n n 

At the lowest order in e we have the system 

A ^-(l + Ao + f) y«°>(x) = 
^% W (*) + |l4 0) (*) + (Ao-l + 2)y(°)(»)=0 



(C.5) 



in which we obviously recognize the Dirac equation for j — 0. The solution of the system have 
therefore the well known expressions, |19( . 

yf ] (x) - iVo v/lTA^ e"$ £vCT-i ( i*i (a, 6, + A & (a + 1, b, £) ) 

(x) = -No V-Ao + 1 e "f f/^- 1 ( iFi(o, 6, - A xF x {a + 1, b, £) ) (C.6) 
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where 



-, , 2 , Vl - a 2 Vl - A 2 - A 
( = 2\/l-Ao i, A = 



a + \/l — A 2 

(C.7) 



while the parameters of the confluent hypergeometric functions are 

A a 



fl = y/l-Q 2 - - 6= 1 + 2^1 -a 2 . (C.8) 

V 7 ! - A 2 

Finally No is the normalizing constant that obviously makes sense only when assuming the spectral 
condition 



^( 1 + (»f? r ) (a9) 

with integer n. Letting 

A = l + a 2 w , (CIO) 

the value for the ground state of the Dirac equation is wq = —.5000066566. 

We now switch the perturbative machine on. In the first place we observe that, according to 
(|C.2f> . the functions </>(°)(x) and 9^ \x) are just obtained by the product with the limiting value 
of the rescaled multiplicative factor, i.e. x~ x (1 + A + a/x) 1 / 2 . In particular, cj)^(x) solves the 
second order equation 

£ ^ - ( x2 ~ M 2 o+Q)2 + 1 -, z x a \ Ta 2 ) ^ (*) = o (en) 

ax \ x 4 (x + x Ao + a) x / 

We then go to the linear terms in the expansion parameter e and, after some simple calculations 
we find a second order equation in <j>i(x) involving Ai and the zeroth order quantities. As it should 
be, according to the standard scheme of perturbation theory, the second order self-adjoint operator 
acting on 4>i(x) is the same as the one acting on 4>o(x) in the zeroth order equation l|(J.ll[) . Indeed: 



d 2 m /x 2 -(xA + a) 2 3 a 2 \ , m , , 

ax z \ x z 4 (x + x Aq + ay x z / 



1_ j- (3x + xA + a) 2 _ 23 .( a; + a;Ao + a) _ a .2 A ^ {x + xXa + a ) ^°)(x)+ 



(x + xA +a)^ d (0) 
4x 2 / dx 



( Al ( 2x(x + xA + Q ) + ^ ) - ^3 + - A + a) (-a + 2 x + 2 x A )) *<°> (x) 



(C.12) 



Therefore the correction Ai can be calculated from equation l|C.12(l in the standard way: we take 
for 4>o(x) and 0q(x) the solution obtained by l|C.6|) and the spectral condition, we multiply l|C.12(l 
by 4>o(x) and integrate the result from zero to infinity. Letting then Ai = a 2 w\ and carrying over 
the numerical calculations we have indicated, we obtain the value w\ = .232995 10 -4 , so that the 
correction to the ground state of the Hydrogen results in (1 — pn) w\ = 0.253509 10 -7 . This result 
is in excellent agreement with what is found from the numerical integration, as presented in Section 
E| When applied to Muonium, 0, where we have (1 — p m uonium) = 0.0096261, the correction would 
amount to 0.22428 10" 6 . 
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D The numerical scheme 



We outline the numerical methods used to integrate the even and odd fourth-order differential 
systems and to find the spectrum. For the numerical treatment we have rewritten the mass pa- 
rameters in the systems according to the definitions (|5.1I5.3|) and we have adopted the independent 
variable x — m e r. 

As a first remark, we see that both systems are singular at zero and infinity. Therefore, in order 
to start the numerical integration, the appropriate asymptotic solutions at the two singular points 
must be determined by the appropriate developments. We report here below and separately on 
some further details of these solutions. As a general feature, however, it turns out that among the 
four independent solutions that can be determined at each singular point, two of them are divergent 
and would produce not square- integrable wave functions, i.e. outside of the Hilbert space where 
we want to solve the spectral problem. These solutions will therefore be rejected. 

Once the two acceptable solutions in zero and infinity have been determined, the eigenvalue 
will be obtained by adapting to our problem the well known "double shooting method" . More 
precisely, for each of the acceptable solutions both at zero and at infinity we start the numerical 
integration up to a chosen crossing point x c . In x c we impose the matching conditions for each 
component of the wave function. These give rise to the following linear system: 



K lV f l \x c ) + K 2 yf 2 \x c ) = K 3 y\°°' l \x e ) + K 4 y^ 2) (x c ) 



(D.l) 



where j = 1,4. The first superscripts " " and " oo " evidently refer to the singular point where the 
integration started, while the second subscripts " 1 " and " 2 " denote the solution chosen at each 
boundary. Being linear and homogeneous, the system (|U> . 1|1 has nontrivial solutions in K\, K 2 , K 3 
and K4 provided its determinant is vanishing. The condition 
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(D.2) 



gives therefore the equation that has to be solved in order to find the spectral values A. Needless 
to say that when the system decouples in two separate second order differential equations, no 
additional modification to the usual double shooting method is necessary and each equation has 
been solved independently. 

We now describe the behavior at the singular points of the two systems. 

D.l The behavior at zero 

We shall first consider the behavior of the systems at the origin by expanding the equations of 
each system in a neighborhood of x = 0. Accordingly, the solutions will be assumed of the form 

N 

yi ( x ) = x»J2 vl j)xj > (* = 1 ' 4 ) (°- 3 ) 

The procedure is the standard one. We substitute (|D.3|) into the system (|4.7|l and take the different 
orders in x. The lowest order gives the indicial equations that produce four different values for the 
exponent i>, only two of which are acceptable. In correspondence to the indices, some relations for 
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the yf are established. It also turns out that the indices are equal both for the even and the odd 
systems. We have a first index 



vi = -l + \jl - — +](!+.]) (D.4) 



with corresponding relations 



2/f=l, yf = -^, 2/f = 2/f = (D.5) 



a 



and a second index (only for j > 0) 

^ = -1 + ^ / (1 - — |D.(il 

for which 



(°) (°) O (°) 1 (°) / -i I \ m -n 

2/1=2/3 =°: 2/ 2 = 1 » 2/4 = (1 + ^2) (D.7) 

a 

The corresponding (very cumbersome) series have been calculated formally up to the eight 
order and have been tested by substitution in the differential equations with the help of a symbolic 
computer package. They have also been tested numerically, with physical values of the parameters 
and varying value of the coordinate x. In addition, this procedure has furnished an effective 
criterion for choosing the boundary point near zero: in fact we have assumed as acceptable a value 
of x where the ratio of the differential equation evaluated on the solution over the solution itself 
was at most of the order of 10 -12 . 
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D.2 The behavior at infinity 

In order to determine the solutions in the neighborhood of infinity we consider the following 
asymptotic expansion for the components of the solution must be looked in the form (see |2Zj) 

JV 

yi (x) = e-r*x» £ V^x-i, (i=l,4) (D.8) 
j=o 

Contrary to what occurs in zero, also the indices of the asymptotic expansion depend upon the 
eigenvalue A. At the lowest order, indeed, we find 

7 {l-pf 4 A 2 (l-p) 4 [ ' 

and we shall take the positive square root to insure the convergence. The next order determines 
the exponent v, given by 

1 / a A 8 a 2 a \ 

The different choices of the two coefficients yf^ and determine then two different series and 
therefore two linearly independent asymptotic solutions of the system. 



D.3 The wave function 

Once the spectrum has been found, we can proceed to determine the wave function. If A» is one 
of the spectral values, then, in the system IjD.lfl we can fix one of the four parameters K t to unity 
and solve for the remaining ones. For instance for the Hydrogen atom and the 2p!/ 2 state of the 
Dirac spectroscopy, for which A„ = —0.1250020774, we find 

Kx = 40.26667371 , K 2 = .6304491167 , K 3 = .7071066401 , K A = l. (D.ll) 
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Therefore the complete wave function will result in 



y j {x)=K 1 yf 1 \x) + K 2 yf 2 \x), 



for x < x c 
for x > x c 



(D.12) 



The plot of yi(x) ig given in the following figure: 




FIG. 3 The first component of the wave function of the Coulomb Hydrogen atom for the 111 + j = i 1 state. 

E The ground state wave functions of Parapositronium and 
Orthopositronium 

In this last appendix we give some details on the two wave functions of the ground state of the 
Parapositronium and Orthopositronium corresponding to the states io+,j = o, 1 and to —,3 = 1,1, 
degenerate in energy. We finally comment on the calculation of the hyperfine splitting of the 
Positronium ground level. 



E.l The Parapositronium 

The case of the Parapositronium, namely the state with even parity, is comparatively easy, since 
it can be solved by discussing the single second order differential equation l|4.13a|) for the Coulomb 
potential, with p = and j = 0. Its explicit form reads: 



dx 2 



2a 



d 



x ((4 + w a 2 ) x + 2 a) J dx 



((■i + w a 2 ) x + 2 a) 2 
16X 2 



l)yi{x) (E.l) 



Zero and infinity are the two singular points of this equation. At infinity there is a divergent 
and a convergent solution, so that the latter will be assumed. In zero both solutions are divergent, 
and only one of them, the solution we shall accept, is normalizable. The regular solution in zero 
has the following expansion: 



2/1(2;) ^ Ax 

x— >0 



1 



(a 2 + 2- V4- a 2 ) (4 + w a 2 ) 



:(4V4~ 



4) 



0- 



(E.2) 



where A is an arbitrary integration constant. 
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The asymptotic behavior of the solution by taking the dominant term of IjE.lfl . at infinity. The 
solution of the corresponding asymptotic equation reads 



yi (x) „ B exp^ — ^ \f—w y 2 + w x ^ (E.3) 

X — >oc 

where B is an arbitrary integration constant. However for calculating a better approximation to 
the asymptotic boundary conditions we can take the expansion of i|E.l|) up to terms of order x~ 2 . 
In fact the resulting equation can be exactly solved in terms of Whittaker M and W functions. 
Up to the usual integration constant we can fix equal to unity, the solution with regular behavior 
at infinity is: 



W (4 + w a 2 ) x + 2 a / 4 + wa 2 \J\ — a 2 a y/—w y/w a 2 + 8 \ 
Vl(x) ~ m W[— ——===, , x J . E.4 

Once the boundary behaviors are established, the wave function component yi{x) and its first 
derivative y[(x) are computed by standard numerical integration. Finally, in order to determine the 
complete sixteen component wave function, we only have to calculate the eight coefficient functions 
a,i(x), bi(x), Ci{x) and di(x), i = 0, 1, in terms of y\{x) and its first derivative y[{x). Recalling that 
2/2(2;) = 0, and that from 14.7fl and (|4.8|l we have 



/ \ Ixy'Ax) ,„ . 

= z u l 2^ (E-5) 

a + x (4 + w a z ) 
simple substitutions lead to the following expressions: 

/ 1 2x \ fl 2 x \ 

a o(x) = ( o + —77-, — 2 \ 1 o y^ x ) ' ai ( x ) = ( o — 7i~, — 2\ 1 o y^ x ) ' 

V2 x (4 + wa z ) + 2 a J \2 xyi + wa 2 ) + 2a/ 

b (x) = bx(x) = co(x) = ci(x) = 0, d (x) = dx(x) = - y 3 (x) . (E.6) 
The components yi(x) and y?,{x) are respectively represented in the first plot of Fig. 4. 



E.2 The Orthopositronium 



Let us then consider the wave function for the odd state. The system to be discussed is now given 
by l|4.7l - 14.81 - I4.11fl with p = and j = 1 . Considering the series expansion in the origin, we see 
that we have two regular solutions with indices 



1 . 



(E.7) 



1, v 2 = \\2 

4 V 4 

On the other side, the asymptotic expansion gives then solutions presenting the exponential de- 
crease 



exp 



(-1 



or 



'—w \l 2 + w — x 



(E, 



exactly the same as for Parapositronium. 

By applying the matching procedures described in Appendix ID. 31 we find the four components 
yi(x). They given in the second plot Fig. 4. 



Letting S be the matrix 



8 = 



2V6 



(4 + a 2 w)x + 2a 



-V2 


1 


-V2 


-1 


-1 


-V2 


-1 


V2 



(8 + a 2 w)x + 2a 
4^2 

a 

= (2 + aw)x 

4a/2 

3 + a 2 w)x + 2a 



— (2 + QIO I 

4 



+ a 2 w)x + 2a \ 



— (2 + aw)x 
4 

1 + a 2 w)x + 2a 
(2 + aw)x 



(E.9) 
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the explicit relationships among the four components yk(x) and the eight coefficients of the com- 
plete state vector, a,(x), bi(x), Ci(x), di(x), i = 0, 1 are the following: 

a Q {x) = a x {x) = -y x (x) , b (x) = -b x (x) = - y 2 (x) (E.10) 

and 

t (c (x),c 1 (x),d (x),dx(x)) = S t (y 1 (x),y 2 (x),y 3 (x),y 4 (x)) . (E.ll) 




FIG. 4 Left. The components of the ground state of Parapositronium: y\(x) (upper) and — 200^3(2) (lower). 

Right. The four components of the ground state wave function of Orthopositronium: from top to bottom 
-wWi -Usix), 100 y 2 (x) and -100 y x (x). 



E.3 The non relativistic limit 

In order to have a better insight into the solutions of the Para- and Orthopositronium ground 
states, we find it useful to look at the Schrodinger limit of the relativistic system and to show how 
the usual equations for the Coulomb problem are recovered. 

In addition to the obvious substitutions /1 = and M — 2m e , the non relativistic limit on the 
system lj4.7N4.8t needs the use of the "atomic parameters" . This means that we have to adopt the 
following rescaling for the independent variable and for the eigenvalue: 

r=— , A = 2m e (l + — ). (E.12) 



As already noted in 14. Ill , when dealing with the even problem and assuming the coefficient 
functions Ij4.11jl . the system 14.7N4.8I) decouples in two separate subsystems, the first one for yi(z) 
and 2/3(2), the second for y 2 {z) and 1/4(2). Thus, isolating 2/3(2), 

fe(*) = 7 T 2 T T 2 (E.13) 

4 z + z cr w + a z dz 

substituting it into the equation for 2/1(2) and taking the limit a — > 0, we find 

^ Sl(z)+ || ! , l(2)+ ( 2m+ |_i(i+i)) !(l(z)= „, (E . 14) 
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in which we recognize the usual £ = j radial Coulomb equation, that gives the well known expo- 
nential solution yi(z) = exp(— z), when j = and w = —1/2. The analogous treatment for the 
second subsystem leads exactly to the same result. 

For the odd case, assuming the coefficients l|4.14jl . the system resulting from (|4.7N4.8|I does not 
decouple: as a consequence the treatment is just a little bit technically more complicated. In order 
to determine the non relativistic limit we have to look for the equations satisfied by the "large" 
components of the wave function. As shown in Fig. 4, they are yz{z) and yi(z). Therefore, after 
changing to the atomic variables, we isolate yi{z) and 2/2(2) from the last two equations of the 
system and we substitute the results in the first two ones. Taking the limit for a — ► 0, we get the 
following system of two coupled second order differential equations: 



dz 2 



2/3 0) 



2 d 

z dz 



2/3(2) + ( 



2w 



j(j + 1) + 2 



2 d 



~T~2 2/4(2) + - — yi{z) + I 2w + - 



dz 



z dz 



3(3 + 1) 



) vM + vW+i) 2/4(2) = 



) 2/4(2) + \ ^fJU+T) y 3 (z) = (E.15) 



The differential operators of both equations of (|E.15|I coincide. The matrix of the non-differentiated 
terms 



2w 



2 j(j + 1) 

z z 2 



2 VJU + 1) 



is diagonalized by the constant matrix 



2 

2w+ - 

z 



3(3 + 1) 



(E.16) 



T 



VTTT V7 



(E.17) 



-VJ VTTT 

Therefore, letting \ipi(z), ^2(2)) = (T^ 1 ) \yi(z),y2(z)) we obtain the two decoupled equations 



d 2 / \ 2d . , / 
— Vl (z) + ~ — y 1 {z)+ 2 
dz z z dz \ 

d 2 .. 2d . . 
dz- 2lP2{z) + -zd-z V2{z) 



2 (j + 1) (j + 2) 
w H 7. 



z 



2 j (j - 1) 
2w+~- JKJ a ; 



.2 



.3^ 



)^(z) = 0. 



(E.18) 



corresponding to the radial Coulomb equations with £ = j + 1 and t = j — l. 



E.4 Some considerations on the hyperfine shifts 

Comparing the limiting classical wave functions with the results reported in Fig. 4, it appears that 
the relativistic corrections become weaker and weaker for large z, as expected by obvious physical 
reasons. Near zero the bare Coulomb potential is no more sufficient for an accurate description 
of the system and, as observed also in case of the Lamb shift, the radiative corrections due to 
the self energy and to the vacuum polarization should be accounted for. The length that, on a 
physical basis, can be assumed as a scale for the relevance of the mentioned relativistic effects is 
the Compton wavelength of the electron A : it is interesting to remark that the sharp maximum 
of all the components of the Orthopositronium occurs approximately at (v2/2) A e . 

These considerations can be brought to bear to a possible evaluation of the hyperfine shift 
of the ground states of the Para- and Orthopositronium. Indeed the usual and very successful 
derivation of the shifts of the two ground levels is done by a perturbative method in a non relativistic 
approximation. Of course, in view of the discussion of the previous subsection, we are able to obtain 
exactly the same results if we adopt the non relativistic approximation and the same expressions for 
the spin-spin interaction and for the annihilation energy term. The latter are derived as effective 
potentials from a low momentum expansion of the one-photon exchange Feynman diagram, |2M| . 
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and the value of the hyperfine shifts of the ground levels are produced by the sum of two different 
contributions. The first one is the average of 47r/io(7S 2 /3 — 2) 5(r) and therefore implies the 
evaluation of the classical Coulomb wave function at zero, as in the original Fermi calculation 
|28| . This approach cannot be directly used with the wave function we have produced, since the 
relativistic effects naturally lead to a drastic modification of their behavior just in the origin. 
The averaging over the angles of the other term, 6/1q ((Sr)(Sr)/r 5 — S 2 /3r 3 ) , giving a rigorously 
vanishing value in the non relativistic limit, is now negligible only for for large values of r as 
expected, but it is divergent in the origin, even after the radial measure r 2 has been accounted 
for. One could try to determine a phenomenological value of the radial coordinate where to put 
a kind of cut-off both for the evaluation of the wave functions and for the integrations as well. 
However, in order to use the relativistic quantum mechanics for calculating the hyperfine splitting, 
we believe it would be much more appealing to deduce an effective interaction potential without 
assuming a low momentum approximation and including radiative corrections. These ideas will be 
developed elsewhere. 
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